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ABSTRACT 

Two  closely  related  algorithms  are  presented  for 
extrapolating  to  the  limit  of  a  scalar  sequence.  One,  the 
even-epsllon  algorithm,  Is  due  to  Wynn{  It  permits  system¬ 
atic  calculation  of  the  array  of  Shank's  transforms  or, 
equivalently,  of  the  related  Pad!  Table.  The  other,  the 
even-rho  algorithm,  is  closely  related  to  the  first  and  is 
also  based  on  Wynn's  work}  however,  it  has  different 
properties  and  has  not  enjoyed  the  same  theoretical 
development.  Singular  rules  and  near-singular  rules  are 
developed  for  both  algorithms  to  handle  situations  in 
which  adjacent  tabular  entries  are  equal  or  nearly  equal, 
leading  to  zero  or  very  small  divisors.  Computer  pro¬ 
grams  implementing  these  algorithms  are  given  along  with 
sample  output.  An  appreciable  amount  or  historical 
background  material  is  included. 

ADMINISTRATIVE  INFORMATION 

This  research  was  carried  out  under  the  Mathematical  Sciences  Research  Program, 
Task  Area  SR-0140301  (Math  Sciences),  supported  by  the  Naval  Sea  Systems  Command. 


1.  INTRODUCTION 


1.1  PURPOSE  AND  SCOPE 


The  purpose  of  this  report  is  to  present  two  computational  algorithms,  called 
the  even-rho  and  the  even-epsilon  algorithms,  for  accelerating  convergence  of  numer¬ 
ical  sequences  and,  in  addition,  to  present  singular  rules  and  near  singular  rules 
for  handling  cases  where  zero  or  near  zero  divisors  appear  in  the  computations.  Nu¬ 
merical  sequences  arise  in  a  wide  variety  of  contexts  in  applied  mathematics  and 
engineering;  examples  are  the  successive  partial  sums  of  an  infinite  series  and  the 

successive  outputs  of  an  iterative  process  s  ,  =  f (s  ).  Singular  cases  require  spe- 

n+i  n 

cial  treatment  to  avoid  excessive  loss  of  significance.  The  computational  algorithms 
as  given  here  are  new  but  their  derivation  depends  heavily  on  the  work  of  Peter  Wynn. 
The  singular  rules  and  near  singular  rules  appear  here  for  the  first  time. 


.1 


The  even-epsilon  algorithm  is  a  computational  adaptation  of  a  formula  given  by 
Wynn  (1966)**  relating  adjacent  members  of  the  Pad6  table;  this  table  is  a  rec¬ 
tangular  array  of  rational  functions  the  series  expansions  of  which  agree  in  a 
specified  manner  with  a  given  series. 

The  even-rho  algorithm  is  based  on  developments  closely  parallel  to  the  above 

2 

but  stemming  from  the  reciprocal  differences  introduced  by  Thiele  (1909)  to 
facilitate  interpolation  by  means  of  rational  functions.  This  algorithm  makes  its 
debut  here. 

These  two  algorithms,  although  very  similar*. -have  different  properties  and 

behave  differently  when  applied  to  any  particular  numerical  sequence.  Five  examples 

of  numerical  sequences  and  the  results  of  applying  the  even-rho  and  even-epsilon 

algorithms  to  them  are  given  in  the  appendixes.  Computer  programs  are  also  given. 

An  effort  has  been  made  to  include  sufficient  historical  background  to  provide 

perspective  for  the  even-rho  and  even-epsilon  algorithms,  i.e.,  to  show  where  they 

fit  into  the  scheme  of  efforts  to  accelerate  convergence  of  numerical  sequences.  To 

a  limited  extent  this  historical  material  also  provides  background  for  a  companion 

3 

report*  Eddy  (1980)  on  a  very  successful  method  of  accelerating  convergence  of  a 
vector  sequence  generated  in  the  iterative  solution  of  a  system  of  linear  algebraic 
equations. 

1.2  SEQUENCES  AND  SERIES 

Numerical  sequences  and  series  which  must  be  evaluated  numerically  arise  in 
many  contexts  in  applied  mathematics.  A  sequence  may  be  the  successive  outputs  of 
an  iterative  process 

s  -  f(s  ,)  (1.2-1) 

n  n-i 

or  it  may  represent  the  successive  partial  sums  of  a  power  series 

n 

s  "  y  a.t^  (1.2-2) 

n  ^  J 

j-o 

*A  complete  listing  of  references  is  given  on  page  53« 
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Conversely,  with  any  numerical  sequence  Sq,  s^»  82*...  there  can  be  associated  a 
power  series 


S<t)  »  Bn  +  /,  {8j-81_,)tj 

J-l 

having  partial  sums 


Sn(t)  "  80  +^7(8j“8j-l)tJ 

which,  for  t  ■  1,  collapse  to  the  original  sequence: 

S  (1)  *  s 
n  n 


(1.2-3) 


(1,2-4) 


(1.2-5) 


Such  a  sequence  or  series  may  converge  with  sufficient  rapidity  to  be  computa¬ 
tionally  useful  as  it  stands;  it  may  converge  so  slowly  that  acceleration  techniques 
are  required;  or  it  may  diverge  so  that  special  techniques  are  required  to  determine 
the  antilimit  from  which  it  is  diverging. 


1.3  CLASSICAL  BACKGROUND 

Ever  since  infinite  series  came  into  use  about  the  time  of  Newton,  the  problem 
of  accelerating  the  convergence  of  slowly  converging  series,  or  of  assigning  a  value 
to  a  divergent  series,  had  to  be  dealt  with.  Various  methods  for  doing  so,  deve¬ 
loped  during  the  18th  and  19th  centuries,  are  discussed  in,  for  instance,  Knopp's 

4 

classic  treatise,  Knopp  (1951),  especially  chapters  8,  13,  and  14,  and  also  in 
Kline  (1972)  ,  chapters  20  and  47.  Among  the  best  known  of  these  methods  is  the 
Euler  transformation 


00 

y* (-l)n  a 

n 

n"0 


(4na  )/2n+1 


o 


(1.3-1) 
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nth 

wher n  4  is  the  n  forward  difference  of  the  e's  beginning  at  a^*  Another  well- 

known  method  which  is  still  in  use  is  the  Euler-Maclaurin  summation  formula 

JT'f(n)  -  J  f(x)dx-{  (f(N)-f(O))  4^  (f  (N) (2N'1  }-f  (0)  (2N“l } )  (1.3-2) 

n"0  0  n"l 

where  the  are  the  Bernoulli  numbers. 

Toward  the  end  of  the  I9th  century  two  new  methods  of  dealing  with  the  conver¬ 
gence  problem  were  developed.  The  first,  called  summability  theory  and  associated 
primarily  with  the  names  of  CesAro  and  Holder,  employed  various  "means"  v.'hich  were 
linear  combinations  of  successive  partial  sums.  These  means  proved  to  be  of  great 
theoretical  value  but  seem  to  be  not  very  useful  computationally.  The  second  new 
method  involved  the  use  of  rational  functions,  either  as  continued  fractions  or  as 
the  ratio  of  two  polynomials}  it  has  largely  supplanted  other  methods  for  acceler¬ 
ating  convergence. 

The  way  was  opened  by  Stieltjes  who,  in  a  series  of  papers  beginning  in  1889, 
exploited  the  idea  of  converting  the  infinite  tail  of  a  series  (the  part  left  after 
removing  a  partial  sum)  into  a  continued  fraction.  Since  the  domain  of  convergence 
of  a  continued  fraction  is  quite  different  iron  that  of  a  power  series,  this  tech¬ 
nique  often  yields  useful  reults.  See  Wall  (1948)^. 

Soon  thereafter  came  PadA  who,  in  his  1892  thesis,  studied  the  relations 

between  a  given  power  series  and  the  rectangular  array  of  rational  functions  (later 

known  as  the  PadA  table)  related  to  it  as  follows:  the  rational  function  in  the 
t  Vi  tVi 

p  row  and  q  column  has  numerator  of  degree  p  and  denominator  of  degree  q  and 
its  power  series  expansion  agrees  with  the  given  series  through  the  term  of 
degree  p  4  q.  (Actually,  the  array  treated  by  PadA  was  the  transpose  of  this  one, 
but  this  description  conforms  to  current  usage.) 

To  evaluate  the  limit  of  convergent  series,  or  the  antilimit  of  a  divergent 
series,  one  evaluates  the  array  of  functions  in  the  PadA  table  and  looks  for  conver¬ 
gence  down  the  columns,  q  «  constant,  or  along  the  main  diagonal,  p  ■  q. 

For  purely  numerical  work,  the  most  effective  way  of  accomplishing  this  eval¬ 
uation  is  with  the  epsilon  algorithm  (or  the  even-epsilon  algorithm)  described  in 
Section  3. 


Considerable  theoretical  work  has  been  done  on  Pad6  approximation  during  the 
past  two  decades;  the  present  state  of  the  art  is  summarised  in  Baker  (1975).^ 
Earlier  work  is  described  in  the  standard  works  on  continued  fractions;  Perron 
(1929)8  and  Wall  (1948). 6 

By  far  the  best  known  and  most  widely  used  example  of  summation  by  approxima¬ 
tion  by  a  rational  function  is  the  simplest  possible  case:  the  rational  function  is 
merely  the  usual  expression  for  the  sum  of  a  geometric  series  with  the  common  ratio 

determined  from  three  successive  partial  sums.  This  simple  function  has  been  dis- 

9  10 

covered  and  used  by  various  authors,  most  notably  Aitken  (1926)  and  Shanks  (1949) 
and  (1955).**  Because  of  the  second  difference  in  the  denominator,  Aitken  called  it 
the  delta  squared  process.  This  extrapolation  process  can  be  written  in  any  of  the 
following  equivalent  forms: 


s  ,  -  (s  -s  , )2/(s  ,-2s  +s  . ) 

n-1  n  n-l  n+1  n  n-1 


s  -  (s  .,-s  )(s  -s  ,  )/(s  ,,-2s  +s  ,) 
n  n+1  n  n  n-1  n+1  n  n-1 


8  .i  “  (s  ,i“S  )  /(s  ,,-2s  +s  .) 
n+1  n+1  n  n+1  n  n-1 


Sn+1 

sn+l 


n-1 


-  2s 


1 


A  - 
n 


s 

n 


s 

n 


s 

n 


(1.3-3) 


(1.3-4) 


Expression  (1.2-3)  is  the  form  usually  given;  (J.2-4)  is  displayed  because  it  is  a 
special  case  of  the  even-epsilon  algorithm  to  be  described  in  Section  3. 

For  a  comprehensive  but  concise  account  of  acceleration  methods  currently  con¬ 
sidered  to  be  of  interest,  especially  the  epsilon  algorithm,  see  the  lecture  notes 
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of  Claude  Brezinski  (1977)  and  also  his  textbook,  Brezinski  (1978). 


2.  THE  RHO  AND  EVEN-RHO  ALGORITHMS 

2.1  THE  RHO  ALGORITHM 

U 

The  rho  algorithm  waa  so  named  by  Peter  Wynn  who,  in  Wynn  (1956),  the  first 
of  two  important  papers  published  a  few  months  apart,  called  attention  to  the  useful 
ness  of  reciprocal  differences  for  summary  series  and  for  extrapolating  to  the  limit 
of  sequences. 

Reciprocal  differences,  denoted  by  p,  had  been  invented  a  half  century  earlier 
2 

by  Thiele  (1909)  to  facilitate  interpolation  by  rational  functions.  The  inter¬ 
polating  function  appears  as  a  continued  fraction  which  is  the  analog  of  Newton's 
polynomial  interpolation  formula  based  on  his  divided  differences.  Analytical  prop¬ 
erties  (.«.  hese  reciprocal  differences  were  investigated  in  the  ensuing  years  by 
Norlund  whose  famous  text,  Norlund  (1924),*'’  remains  the  definitive  treatment  of  the 
subject.  A  similar  treatment  is  found  in  Milne-Thomson  (1933).^ 

The  reciprocal  differences  corresponding  to  a  given  set  of  arguments  {X^}  and 
function  values  {s  >  are  defined  recursively: 


. ,  X  -  X 
n  n+1  ,  n+m  n 

pm  "  pm-2  n+1  n 

pm-l  pm-l 


(2.1-1) 


n  n  _ 

P„  "  s„»  p-i  s  0 

o  n  -1 


(2.1-2) 


They  are  customarily  displayed  in  the  rho  array  a-  shown  in  Figure  2.1. 

,1 


'i 


P  *>  s, 
0  ’ 


^-p 


P0^S2 


P3=S 
"o  s3 


P*  “  s4 
0  4 


P6  =S 
"  0  s5 


•-P 


,  _X2  X1 
1  '  S2  S1 
2  x3  '  x2  _  - - 

1  “  V  S2 


1  2 
=  P0  + 


X3  '  X1 

P2  •  P1 
r\ 


X4  X2 


3  X4  ’  X3 

P  - - 

1  S4  •  S3 


.  xc  ■  X„ 

P4  =  _L_J 

1  V  S4 


^  =  Pp  +  P3  -  p2 

1  1 


;3  =  p4  + 


P2  =  P0 


X5  ‘  X3 


Figure  2.1  -  The  Rho  Array  (Old  Notation) 
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It  will  be  seen  that  any  particular  p^  depends  explicitly  only  on  three  pre¬ 
viously  calculated  rho's  lying  at  the  other  three  vertices  of  a  rhombus.  However, 
implicitly  pP  depends  on  all  m+1,  of  the  pairs  (X^,  s^),  sn+i ) • • • ^n+m* 


Vh.>- 

Systematic  computation  of  the  rho  array  can  proceed  according  to  either  of  two 
patterns: 


columnwise  -  all  of  the  input  sequence  p 
n 

p^,  etc.,  or 

by  upsloping  diagonals  -  as  each  new  pn  « 
n-1  1 


V 


then  all  of  the  first  column 


s  is  obtained,  calculate  in  turn 
n 


P1  ’  Pn  1*  latter  pattern  has  two  distinct  advantages, 

•  it  reduces  storage  requirements  since  calculation  along  any  upsweeping  diag¬ 
onal  requires  input  from  only  itself  and  the  immediately  preceding  upsweep¬ 
ing  diagonal,  and 

•  convergence  can  be  monitored  in  all  columns  involved  on  each  sweep  so  that 
the  basic  input  interation  can  be  halted  as  soon  as  possible. 

The  property  of  reciprocal  differences  which  makes  them  so  useful  for  ..extrapoK,  ■ 
lating  to  the  limit  of  a  sequence  (or,  equivalently,  for  summing  a  series)  is  that, 
if  the  general  term  of  a  sequence  is  given  by  a  rotional  function  of  the  term  index 
(ordinarily  X^=n)  in  t^e  form 


;  -4 


,  5 


then 


s  = 


a,  Xk  +  a.  ,  Xk_1  +  . . .  +  a,  X  +  a 
k  n  k-1  n  1  n  0 


Xk  +  b.  ,  Xk_1  +  ...  +  b.  X  +  b_ 
n  k-1  n  l,n  0 


(2.1-3 ) 


-.1 


p2k  "  ak 


for  n  =  1 ,  2,  3,  . . 


(2.1-4) 


Since  obviously 


lim  s 

r 

x  -*•<= 
n 


(2. 1-5) 


it  is  clear  that  calculation  of  the  even  ordered  reciprocal  differences  provides  a 
systematic  way  of  finding  the  limit  or  antilimit  of  such  a  sequence.  If  relation 


(2.1-3)  is  not  exact,  but  is  an  increasingly  good  approximation  as  n  gets  larger, 

then  there  will  be  convergence  in  the  2k  column  of  the  rho  array. 

The  rho  algorithm,  in  contrast  to  the  epsilon  algorithm  discussed  in  Section  3, 

seems  to  have  received  very  little  attention  in  the  literature  since  the  publication 
14 

of  Wynn  (1956)  and  convergence  theorems  for  it  are  lacking.  It  has,  however,  been 
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touched  upon  several  times  by  Brezinski  (1971),  (1972),  and  especially  (1977). 

2.2  THE  EVEN-RHO  ALGORITHM 

Since  only  the  even-ordered  reciprocal  differences  are  useful  for  the  summation/ 
extrapolation  process,  it  would  be  desirable  to  develop  a  recursion  relation  involv¬ 
ing  them  alone.  This  Indeed  can  be  done  and  the  result  is  Equation  (2.2-4)  which,  to 
the  best  of  the  author's  knowledge,  appears  here  for  the  first  time  in  print.  The 
novelty  is,  however,  trivial  since  the  pattern  of  elimination  which  leads  to  this 
recursion  relation  is  precisely  the  one  employed  by  Wynn  (1966)  to  develop  the  even- 
epsllon  algorithm  from  his  earlier  epsilon  algorithm.  In  order  to  make  this  report 
more  .self  contained,  the  derivation  is  given  here  for  the  even-rho  algorithm. 

Consider  a  rhombus  shaped  portion  of  the  rho  array  as  shown  in  Figure  2.2. 


Figure  2.2  -  A  Portion  of  the  Rho  Array 


The  rho  recursion  relation  (2.1-1),  now  for  convenience  rewritten 


A*  •  *  . 

p2k-l  + 


X2k+n+l  ~  Xn 


(2.2-1) 


gives 


n  n+1  /v  v  W/  n+1  n 

p2k+l  "  p2k-l  "  (X2k+n+l"V/(p2k  -P2k) 
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,  n+l  n  v  X2k+n+l  ^2k+n  \-l 

P2k-1  P2k-1'  n+l  n  n  n-1 

p2k  "  P2k  P2k  P2k 

/■  . 

Now  Equation  (2.2-1)  can  be  used  to  replace  the  ordinary  difference  on  ti.e  left 
with  reciprocal  differences: 

*2k+n+l  "  Xn-1  _  X2k+n  "  Xn  =  X2k+n+l  "  Xn  _  *2k-Hi  ~  Xn-1  „ 

n-1  n  n  n+l  n+l  n  n  n-1  '  '  ' 

p2k+2  “  p2k  p2k  "  p2k-2  p2k  “  p2k  p2k  “  p2k 

Note  that  on  the  left  side  the  differences  lie  in  a  horizontal  direction*  whereas  on 
the  right  side  they  lie  in  a  vertical  direction. 

For  ordinary  use  in  summation  of  series  or  extrapolation  of  sequences  it  is 
natural  to  take 


,  n  n-1  . 
p2k+l-p2k+l 


i 


r  ; 


1 


>  ) 


Xn  “  n  (n  »  1,  2,  3,  ...) 


(2.2-3) 


This  substitution  yields  on  the  even-rho  algorithm 


2k+2 _ 2k  _  2k+l  _  2k+l 

n-1  n  n  n+l  n+l  n  n  n-1 
p2k+2  “  p2k  p2k  “  p2k-2  p2k  "  p2k  p2k  “  p2k 

pn  -  sn  (n  “  1,  2,  3,  ...)  (k  *  0,  1,  2,  3,  ...) 


It  is  not  necessary  to  specify  the  since  the  only  term  in  which  they  would  ; 

appear  has  k*0  in  the  numerator.  (  \ 


There  is  now  a  cross  pattern,  shown  in  Figure  2.3,  which  corresponds  to  the 
rhombus  pattern  associated  with  the  ordinary  rho  algorithm. 


jn  1 
2k 


>"  +  1  - 
2(k  •  1) 


2k 


>"  +  1 
2k 


I"’1 
2 (k  +  1) 


Figure  2.3  -  The  Even-Rho  Pattern  (Old  Notation) 


2.3  NEW  NOTATION 

The  notation  p^.  used  so  far  has  corresponded  to  that  used  by  Wynn  in  his 
later  papers  ori  the  epsilon  algorithm  and  also  by  Brezinski.  The  lower  index,  m  or 
2k,  denptes  the  column,  and  the  upper  index,  n,  denotes  the  downsloping  diagonal  on 
which  it  lies.  The  latter  thus  emphasizes  the  first  member  of  the  input  sequence, 
s  ,  which  enters  into  the  calculation  of  p"  (see  Figure  2.1). 

T1 

However,  for  a  computational  algorithm  which  sweeps  through  the  array  along 
upsloping  diagonals,  it  is  much  more  convenient  to  have  the  upper  index  denote  the 
ups loping  diagonal.  Hence  it  is  constant  along  each  sweep,  and  thus  emphasizes  the 
last  member  of  the  input  sequence,  s^,  which  enters  into  the  current  calculation. 
Henceforth  the  upper  index  denoting  the  upsweeping  diagonal  will  be  p  instead  of  n; 
they  are  related  by 


p**n  +  2k  (2.3-1) 

In  this  notation  the  even-rho  array  appears  as  in  Figure  2.4. 

The  even-rho  recursion,  relating  the  five  elements  lying  on  a  cross  pattern, 
appears  as  in  Figure  2.5. 
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Equation  (2.2-4)  becomes 


-  2,(k.tlJ _ 

p+1  p 

p2(k+l)  "  °2V 


2k 

p  p-1 

p2k  "  p2(k-l) 


2k+l 

,P+1  -  0P 


p2k  '  p2k 


(2.3-2) 


(p  ■  2,  3.  4,  ...»  k  “  0,  1.  2,  ...  [ (p— 1)/2J) 

Pq  "  Sp  ^P  "  3,  ...) 

The  new  notation  has  an  additional  advantage — mainly  esthetic,  to  be  sure — 
in  that  the  upper  indices  in  Equation  (2.3-2)  behave  properly  instead  of  running 
backward  in  the  horizontal  differences  a3  they  did  in  the  old  notation  of  Equation 
(2.2-4). 

2.4  COMPUTATIONAL  ALGORITHM 

The  two  terms  on  the  left  side  of  Equation  (2.3-2)  have  the  same  structure, 
that  of  principal  parts  of  reciprocal  differences  taken  in  a  horizontal  direction. 
Likewise  the  two  terms  on  the  right  are  principal  parts  of  reciprocal  differences 
taken  in  a  vertical  direction.  Each  such  term  is  uBed  twice,  once  on  each  of  two 
successive  sweeps  along  upsloping  diagonals.  It  is  therefore  advantageous  to  define 
the  auxiliary  quantities 


P  p-1 

if  -  of 


(2.4-1) 


2 (k— 1 > 


(2.4-2) 


p2k  -  p2k 


In  terms  of  these  quantities,  the  equations  for  calculating  the  along  an 
upsloping  diagonal  indexed  by  p,  as  in  Figure  2.4,  become 


12 


yP  -  (P  -  2,  3,  5,  ...) 


P  P“J 
»2k  '  "2k 


“mi  ■  Hr‘ +  vk  -  vr‘  <»  ■  j> *• »•  •••> 


(2.4-3) 


(2.4-4) 


’2wu  ’  “K1  :p  ■  3- *■  s.  •••)  «•«> 


(K  -  0,  1,  2,  ...  [(p-D/2),  p  -  2,  3,  4,  ...) 


(2.4-6) 


with  initial  conditions 


pP  -  s 
o  p 


Hp  -  0 
o 


(p*  1 »  2,  3,  ...) 


(2.4-7) 


Since  calculations  along  any  upsloping  diagonal,  indexed  by  p,  require  as  input 
only  previously  computed  values  for  the  same  p  and  also  for  p-1,  it  is  necessary  to 
maintain  six  arrays:  three  for  the  currently  computed  values  of  V,  H,  and  rho,  and 
three  similar  arrays  from  the  previc  -s  (p-1)  sweep.  Their  maximum  length,  L,  is  re¬ 
lated  to  P,  the  maximum  value  of  p,  by 


L  -  [(P+D/2] 


where  [x]  means  the  greatest  integer,  n,  satisfying  n  <_  x  <  n+1. 


3.  THE  EPSILON  AND  EVEN-EPSILON  ALGORITHMS 

3.1  SHANKS'  TRANSFORMS  AND  WYNN'S  EPSILON  ALGORITHM 

The  epsilon  algorithm  was  so  named  by  Peter  Wynn  who  introduced  it  in  Wynn 
19 

(1956)  as  a  practical  method  of  calculating  the  array  of  transforms  e  (S  ) 
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discussed  by  Shanks  in  his  thesis,  Shanks  (1955).  The  latter  is  credited  by  Gragg 
20 

(1972),  page  2,  as  being  one  of  two  principal  Bources  of  stimulus  for  modern  inter¬ 
est  in  Pad4  approximation  and  related  topics.  Actually,  an  earlier  version,  Shanks 
(1949),  had  attained  unusually  wide  circulation  £of  an  internal  memorandum  and  was 
widely  known  among  numerical  analysts  in  the  early  1950's. 

Shank's  considered  sequences  of  the  form 

m 

■.  •  «.  'Z  Vi; 

k-1 

which  arose,  for  example,  in  studying  the  decay  of  a  mixture  of  radioactive  sub¬ 
stances.  If  all  r^  are  distinct  and  satisfy 

|rkl  <1  (3.1-2) 


then 


lim  8 


n-»® 


n 


a 

o 


(3.1-3) 


1 

I 

j 


so  that  a  is  the  quantity  of  interest  in  this  context, 
o 

Two  properties  of  the  sequence  (3,1-1)  should  be  noted: 

(1)  each  is  the  sum  of  the  n  partial  sums  of  m  different  geometric  series; 
m 


(2) 


s 

n 


a 

o 


has  the  form  of  the  general  solution  of  a  homogenous 


linear  finite  difference  equation  with  constant  coefficients 


m 

^7Cj<Sn+j“ao)  "  °  (n  "  l>  2>  3»  •••>  (3.1-4) 

j-o 

where  the  r^  are  the  roots  of  the  characteristic  equation 
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for  estimating  aQ,  the  limit  (or  antilimit)  of  the  sequence  in  terms  of  2m+l  succes¬ 
sive  elements,  s  ,  s  . ..  s  .  Here,  n  denotes  the  starting  point  in  the 
n  PTi  n+im 

sequence  and  m  denotes  the  order  of  the  approximation  as  in  Equation  (3.1-1).  In 

particular,  for  ra"l  this  yields  the  well  known  Aitken  delta-squared  extrapolation 

formula  already  shown  in  Equation  (1.3-3), 

The  expressions  (3.1-7)  given  by  Shanks,  the  ratio  of  two  determinants  of  Kankel 

type,  become  practically  useless  for  direct  computation  even  when  in  is  still  a  small 

1  o 

Integer,  e.g.,  m"4.  This  difficulty  was  soon  surmounted  by  Wynn  (1956)  "  who  showed 
that,  with  the  aid  of  some  new  intermediate  quantities,  a  very  simple  recursion  rela¬ 
tion  holds. 

He  set 


C2m  '  W 


*«  .,  ■  1/e  (s  ,-s  ) 
2m+l  m  n+1  n 


(3.1-8) 


and  showed  that  they  satisfy 


en+l  +  - l- - 

k-1  n+1  n 

£k  -  Ek 


k«0,  1,  2,  3, ... 
n«0,  1,  2,  3,... 


(3.1-9) 


en  ■  s 
o  n 


This  is  the  epsilon  algorithm;  it  is  very  similar  to  the  rho  algorithm.  Equation 
(2.1-1).  Correspondingly,  there  is  the  epsilon  array  shown  in  Figure  3.1. 
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Sti. 


Figure  3. 1  -  The  Epsilon  Array  (Old  Notation) 


Systematic  computation  of  this  array  can  proceed  either  columnwise  or  along  up- 

sloping  diagonals  as  was  discussed  for  the  rho  algorithm  in  Section  2.1. 

The  relationship  between  the  array  of  transforms  (3.1-7)  and  the  rational 

functions  of  the  Pad6  table  was  pointed  out  in  Shanks  ( 1 94 9 ) 1 ^  and  ( 1 955) ; 1 ^  the 

same  relationship,  but  in  the  language  of  the  epsilon  array,  was  treated  in  Wynn 
2 1  th 

(1961).  Specifically,  let  s^  be  the  n  partial  sum  of  a  power  series 


(3. 1-10) 


Substitute  Equation  (3.1-10)  into  the  expression  (3.1-7)  for  the  e  (s  )  trans- 

in  n  ^ 

form;  then  in  both  numerator  and  denominator  multiply  the  first  column  by  z  ,  the 

ift**  1  0 

second  column  by  z  ,  ...  the  last  column  by  z  ;  cancel  out  the  powers  of  z 


in  corresponding  rows  of  numerator  and  denominator.  What  is  left  is  very  clearly 
the  ratio  of  a  polynomial  of  degree  (n+m)  to  a  polynomial  of  degree  m: 


3.2  THE  EVEN-EPSILON  ALGORITHM 

As  with  the  rho  algorithm,  since  only  even-ordered  epsilons  are  useful  for 
extrapolation,  it  would  seem  desirable  to  eliminate  the  odd-ordered  epsilons  and 
arrive  at  a  recursion  relation  involving  only  even-ordered  epsilons.  This  elimina¬ 
tion  was  accomplished  by  Wynn  (1966).  The  details  will  not  be  repeated  here 
since  they  have  already  been  given  in  Section  2,2  on  the  rho  algorithm.  Likewise, 
the  arguments  for  changing  to  a  new  notation,  already  given  for  the  rho  algorithm 
in  Section  2.3,  will  not  be  repeated  here.  Suffice  it  to  show  in  Figure  3.2  the 
even-epsilon  array  in  the  new  notation! 
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P  k-0 


1  V‘l 


2  <o-5a 


3  (V‘3 


4  e0*’4 


6  1  0  *  *5 


6  ^ 
0  s6 


7  *0**7 


Figure  3.2  -  The  Even-Epsilon  Array  (New  Notation) 


Any  five  entries  in  this  table,  lying  on  a  cross  pattern  as  indicated  in  Figure 


21k  +  1) 


Figure  3.3  -  The  Even-Epsilon  Pattern  (New  Notation) 


are  related  according  to  Wynn's  formula  (Wynn  (1966),  page  266) 


1 _ 1  _  1  _  1 

p+1  P  "  P  P-1  P+1  P  "  P  P-1 

c2(k+l )  “  e2k  e2k  "  E2(k-1 )  E2k  "  e2k  e2k  "  c2k 


(3.2-1) 
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(p  •  2,  3,  k**0,  1,  2, . . .  ( (p-l)/2J ) 

*o  "  V  '-2  "  *  (P  "  1»  2*  3»***) 


This  is  Che  recursion  relation  among  his  transforms  which  Shanks  came  close  to  but 
never  actually  attained.* 

Comparing  this  expression,  Equation  (3.2-1),  with  that  for  Aitken's  delta- 
squared  process  (Shanks'  e ,  process)  shown  in  Equation  (1.3-4),  we  see  that  a 
quantity  E^  given  by 


l 


(3.2-2) 


is  also  given  by 


1 _ 

p+1  p 

e2(k+l)  "  e2k 


1 

P  p-1 
6 2k  ’  C2(k-1) 


(3.2-3) 


That  is  to  say,  at  a  given  location  in  the  even-epsilon  array  (Pad6  table),  the 
same  result  is  obtained  by  applying  the  delta-squared  (e^)  transform  either 
vertically  or  horizontally. 

3.3  COMPUTATIONAL  ALGORITHM 

The  two  terms  on  the  right  side  of  Equation  (3.2-1)  have  the  same  structure, 
that  of  principal  parts  of  the  next  higher  odd-ordered  epsilcns  (see  the  epsilon 
recursion  relation,  (Equation  (3.1-8)),  the  differences  being  taken  in  a  vertical 


*Oral  communication  from  Shanks. 
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direction.  A  similar  remark  applies  to  the  two  terms  on  the  left  side,  the  differ¬ 
ences  being  taken  in  a  horizontal  direction.  This  observation  leads  to  the  defini¬ 
tions 


P  p-1 
£2k  “  C  2 (k— 1 ) 


(3.3-1) 


-P  _  cP”1 

l2k  e2k 


(3.3-2) 


In  terms  of  these  quantities,  the  equations  for  calculating  the  along 
an  upsloping  diagonal  indexed  by  p,  as  in  Figure  3.2,  become 


P  _  J>“1 
'2k  2k 


2,  3,  4, .  •  • ) 


(3.3-3) 


“T1  +  vk  -  v£_1 


(p  =  3,  4,  5,...) 


'2(k+n  *  e2k‘  +  3~~  (P  -  3,  4,  5,...) 
Hk+1 


(3.3-4) 


(3.3-5) 


(k»0, 1,2,3,... i\p-l)/2],  (p  -  2,  3,  4,...) 


(3.3-6) 


with  initial  conditions 


(p  «  1,  2,  3,...) 


(3.3-7) 


Note  that,  for  k«0,  Equations  (3.3-3)  and  (3.3-4)  together  yield  the  Aitken  delta- 
squared  transform,  Equation  (1.2-4). 


fjrs1 rrrrfrv 


rjjp  ipi'i-; 


As  in  the  case  of  the  even-rho  algorithm,  it  is  necessary  to  maintain  two  arrays 
for  each  of  V,  H,  and  epsilon;  one  each  for  the  current  values  of  V,  H,  and  epsilon 
along  the  upsloping  diagonal  indexed  by  p,  and  one  each  ,fbr  the  last  preceding  diag¬ 
onal  indexed  by  p-1.  Their  maximum  length.  A,  is  again  related  to  P,  the  maximum 
value  of  p,  by 


L  -  [(P+D/2] 


where  [x]  means  the  greatest  integer,  n,  satisfying  r  <_  x  <  n-KL. 


(3.3-8) 


4.  SINGULAR  AND  NEAR-SINGULAR  RULES  .. 

4.1  NEED  FOR  SINGULAR  AND  NEAR-SINGULAR  RULES 

Computational  algorithms,  such  as  the  even-rho  and  even-epsilon  algorithms  and 
also  their  ancestors  the  ordinary  rho  and  epsilon  algorithms,  in  which  division, 
plays  a  major  role,  are  vulnerable  to  troubles  arisipg  from  attempted  division  by 
zero  (the  singular  case)  or  by  a  number  which  is  very  small  in  magnitude  (the  near¬ 
singular  case).  These  situations  are  nuisances  but  need  not  stop  the  calculations; 
instead,  one  can  use  special  formulas  called  singular  rules  for  a  zero  division  or 
near-singular  rules  for  a  near  zero  division.  . 

The  near-singular  rules  will  be  derived  for  the  even-rho  algorithm.  Trivial 
modifications  then  yield  the  near-singular  rules  for  the  even-epsilon  algorithm. 

In  either  case  the  singular  rules  are  then  an  obvious  limiting  case  of  the  near¬ 
singular  rules. 

Both  singular  and  near-singular  rules  for  his  epsilon  algorithm  were  given  by 
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Wynn  (1962)  and  (1963).  He  also  gave  a  related  discussion  for  the  even-epsilon 
algorithm  in  Wynn  (1966).* 


4.2  NEAR-SINGULAR  RULES  FOR  THE  EVEN-RHO  ALGORITHM 

A  glance  at  the  computational  algorithm,  Equations  (2.4-3)  -  (2.4-7),  reveals 
two  places  in  which  a  small  division  may  arise:  the  first  is  in  the  computation  of 
v£,  the  second  in,. the  computation  of  In  both  cases  the  effects  of  near 

singularity  are  felt  at  several  neighboring  points  encountered  later  ("down-stream") 
in  the  calculations.  It  is  interesting  to  note,  in  both  cases,  instances  wherein 


very. large  numbers  cancel  out  by  subtraction;  they  cancel  exactly  to  within  the  accu¬ 
racy  of  the  approximations  used.  These  are  underlined  in  the  formulas  in  which  they 


occur. 


Treatment  is  by  a  first  order  perturbation  method:  the  small  divisor,  d,  is 
assumed  to  have  the  following  properties:  v  ' '' 


(D  Id  I  «  1 

(2)  1 /  |  d  |  >>  any  other  terms  addedto  it 

2  2 

(3)  d  »0,  i.e.,  terras  in  d  are  negligible 


(4.2-1) 


The  simpler  case  will  be  discussed  first.  Suppose  that,  in  the  calculations 
associated  with  the  index  values  p  and  k, 


P  „  uP"1  +  vp  -  vP"1  . 
k+1  k  k  k 


where  d  satisfies  the  criteria  just  given.  Then,  from  Equation  (2.' 4-5), 


(4.2-2) 


P  p-1  2k+2 

p2(k+l)  "  p2k  +  p 

Hk+1 


(4.2-3) 


and  from  Equation  (2.4-3) 


_ 2k+3 _ 

P  P"1 

p2(k+l)  p2 (k+1 ) 


1  2k+2 


(4.2-4) 


On  the  next  following  upsloping  diagonal,  indexed  by  p+1, 


k+1  p+1  p 

p2(k+l)  "  p2(k+i) 


(4.2-5) 


+  vp+1  -  VP  ~ 


k+1  k+1 


d  -  2 


(4.2-6) 


P+1  p  2k+4  _  p-1  2k+2 

p2 (k+2 )  =  p2(k+l)  up+l  ~  p2k  d 

tc+2  - 


(2k+4)d 


-  p^1  +  0  (4.  2-7) 
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These  results  are  summarized  graphically  in  Figure  4.1. 


Figure  4.1  -  First  Case  of  Near-Singular  Rules 
for  the  Even-Uho  Algorithm 


The  corresponding  singular  rules  are  now  trivial:  for  d  =■  0, 


k+1 


p+1 

k+2 


0 


(4.2-8) 


P2 (k+1 )  "  “  (4. 2-9) 

p+1  p-1 

P2 (k+2)  *  p2k  (4.2-10) 
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A  more  complicated  situation  arises  if  a  newly  calculated  p!?,  is  almost  equal 

1  4.JS.  , 


to  p5k  directly  above  it  in  the  even-rho  array: 

••'•v  ,•  (1,  • ;  .... 

s.  •  ,:-;w  ■  .  v:  ,2k,  °2k  d 

V;,  r. v 

'^iH>fesi.';iLiwBisd Lately ••  •* "•  ' 

v  vp-yic  =;j  •?  .',y|  y>“’  •!•■'■'•■  i|  -‘  ■■  •  ny .  . 

:  •  f- ' rr. V,  •  vp  „  2k+i 

Vk  d 


(4.2-11) 


(4.2-12) 


i-.--.V-.  fc  .  ,  d  - 

;;  . (a. .  iv  ’w\2-;  :  ..  ; 

f  ;•  • -  ••  :!  ;  vjy.'  .s  ......  ....  .  .  .. 

and,  under  the  previous1 '^sisumtjlttona  :  of  Eqikktihn  (4.2~1)  about  d,  leads  i.o  the  follow- 

i.  V-  *  !  *  ;  ;  • : \  •'  ■ '.(  :  "  F 


ing  approximations: 


Hk+i  -X" 


pp  i. ...  ii 

P2(k4-1) 
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d 
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■  2k 

\2^1  ) 

(4.2-13) 
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(4. 2-1 4b) 
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2k  P 
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(4.2-Hc) 


Along  the  next  upsloping  diagonal,  .indexed  by  p+1,  the  approximations  become 


Hl+1  “  Hf  +  Vf+l  -  vu  *  “ 

k+i  k  k  k  k  d 


(4.2-15) 


Pp+2  .  pp+1  + 

P2 (k+2)  p2 (k+1 ) 


(4.2-21a) 
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P2(k+1)  ‘  „p+2  +  p+2  "  P2(k+1)  yp+2 

k+1  k+1  vk+l 


(4.2-21b) 


Again,  the  corresponding  singular  rules  are  considerably  simpler:  for  d  ■  0, 


jP"1  .  PP  -  PP  -  PP+1 
2k  P2k  p2 (k+1 )  p2 (k+1) 


■  P.  say; 


(4.2-22) 
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P2 (k+2)  yP 
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2  (k+2) 


P  + 


21k+2) 

vp+r 
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(4.2-24) 


(4.2-25) 


In  Figure  4.2,  the  very  heavy  dots  indicate  those  tabular  entries  of  rho  or 
epsilon  which  are  involved  in  the  second  case  of  the  near-singular  rules.  Heavy 
bars  between  the  heavy  dots  indicate  very  large  values  of  the  corresponding  H  or  V; 
these  values  become  infinite  in  the  singular  case. 


Figure  4.2  -  Second  Case  of  Near-Singular  Rules  for  the 
Even-Rho  and  Even-Epsilon  Algorithms 


4.3  NEAR-SINGULAR  RULES  FOR  THE  EVEN-EPSILON  ALGORITHM 

The  opening  remarks  in  Section  4.2,  and  in  particular  the  restrictions  of  Equa¬ 
tion  (4.2-1)  on  the  small  divisor,  d,  apply  as  well  to  the  even-epsilon  algorithm. 

The  only  difference  is  that  now  H  and  V  are  defined  by  Equations  (3.3-1)  and  (3.3-2), 
respectively,  and  the  whole  algorithm  is  given  by  Equations  (3.3-3)  thru  (3.3-7). 

To  facilitate  comparisons  between  the  corresponding  near-singular  rules  for  the  even- 
rho  and  even-epsilon  algorithms,  corresponding  formulas  bear  corresponding  numbers. 
Equation  (4.3-1)  is  missing  but  it  would  be  identical  to  Equation  (4.2-1). 
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In  the  first  case  of  near  singularity,  let 

\pt,  ■  ‘C1  +  v£  -  r  ■ * 

where  d  meets  the  restrictions  of  Equation  (4.2-1).  Then  Equation  (3.3 
definition  of  H  for  the  even-epoilon  algorithm,  gives 


Ep':' "  „  eP*1  .  _ L.»  I 

*2 (k+1)  2k  +  „p  d 
Hk+1 


and  Equation  (3.3-2),  the  definition  of  V,  gives 


k+1  ,p  p-1 

2 (k+1 )  ~  2(k+l) 


as  d 


On  the  next  upsloping  diagonal,  indexed  by  p+1, 


, _ i _ - _ 

k+1  p+1  p 

e2(k+l)  ”  E2(k+1) 


as  -  d 


«£!  '  «L,  +  *fl!  -  Vj ,  .  *  d  -  2d 

k+2  k+1  k+1  k+1 


-  d 


J _ 


.p+1  P  , 

'2 (k+2)  2 (k+1)  uP+l 

Hk+2 


p-1  ,1  1  p-1 

er  +  —  —  —  =  er 
2k  d  d  2k 


These  results  are  summarized  in  Figure  4.3: 


(4.3-2) 
1 ) ,  the 

(4.3-3) 

(4.3-4) 

(4.3-5) 

(4.3-6) 

(4.3-7) 
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This  configuration  gives  immediately 
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i  vd 


(4. 3-12) 


and,  under  the  stated  assumptions  of  Equation  ( A . 2—1 )  about  d,  leads  to  the  following 
approximations.* 


Cl 


2  (k+1)  2k 


Hu”1  +  v£  "  Vup_1  • 
k  k  k  k  d 


-  C1  +  -  cP-1  +  d  »  eP 
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'2k 


(4.3-13) 


(4.3-14) 


k+1 


Along  the  next  upsloping  diagonal,  indexed  by  p+1,  the  approximations  become 


"&!  - 'f 1  -  --»£*-* 


(4.3-15) 
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P-1 

2k 


VP+1 

k+1 


1 


P+1  P 

‘2  (k+1 )  '  e 2 (k+1 ) 


e2k  “  d  "  C2k 


J. 

d 


H,P+^  =*  HP  +  Vp+*  -  VP  -  VP 
\+2  "k+1  k+1  k+1  \+l 


p+1  p  1  p  _ 

e 2 (k+2 )  "  C2(k+1)  +  Hp+1  "  E 2 (k+1 )  "  „p 


1 


k+2 


k+1 


(4.3-16) 
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(4.3-19) 
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Still  another  upsloping  diagonal  is  also  affected,  that  Indexed  by  p+2i 


UK  '  * 

”k+2 


vp+2 . 

vk+l 


i,P+1 

k+1 


-i 


vP+2 
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+  I 
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"k+2  Vk+1 


The  corresponding  singular  rules,  for  d  *0,  are 


eP”1  .  CP  „  eP  ,  £P+1 
2k  2k  2  (k+1)  2  (k+1) 


e,  say; 
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(4.  3-20) 


(4.3-21) 
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5.  IMPLEMENTATION  OF  THE  ALGORITHMS 


5.1  COMPUTER  PROGRAMS 

In  order  to  carry  out  computational  experiments,  a  small  package  of  computer 
programs  was  developed.  The  even-rho  algorithm,  as  described  in  Section  2.4,  was 
Implemented  as  subroutine  EVRHO  in  CDC  FORTRAN  Extended  (essentially  FORTRAN  IV). 
Similarly,  the  even-epsilon  algorithm  described  in  Section  3.3  was  implemented  as 
subroutine  EVEPS.  In  addition,  a  merged  version  of  these  subroutines,  called 
TANDEM,  was  written  to  take  advantage  of  the  fact  that  they  were  almost  identical, 
differing  only  in  two  lines  of  coding.  The  price  was  to  rename  most  of  the  inter¬ 
mediate  quantities  and  to  double  the  amount  of  storage  required  which  was  trivial. 
These  subroutines  did  not  incorporate  the  singular  rules  or  near-singular  rules 
described  in  Section  4  because  they  were  written  before  the  need  for  near-singular 
rules  was  perceived  and  because  the  necessary  modifications  would  have  been  relatively 
complicated.  Protection  was  built  in,  however,  which  aborted  computations  whenever 
a  divisor  became  too  small. 

Each  sequence  to  be  extrapolated  requires  the  writing  of  a  brief  subroutine  SEQU 
to  calculate  the  successive  numbers  of  the  sequence.  Five  examples  are  listed  in 
Appendix  A.  In  each  case  the  sequence  is  identified  in  the  included  comment  cards. 

In  addition,  there  was  an  executive  routine  TDMCHK  to  perform  the  usual  chores 
of  calling  subroutines  and  handling  output. 

Listings  of  all  of  these  programs  are  given  in  Appendix  A. 


5.2  NUMERICAL  EXAMPLES 

Corresponding  to  each  of  the  sequence  generating  subroutines  SEQU  listed  in 
Appendix  A,  is  an  output  page  in  Appendix  B  which  shows  how  EVRHO  and  EVEPS 
performed  in  that  sequence.  As  might  be  expected  from  the  theoretical  discussions 
in  Sections  2.1  and  3.1,  these  two  algorithms  performed  rather  differently,  one  being 
superior  to  the  other  for  any  particular  sequence  but  neither  being  consistently 
superior  for  all  sequences.  In  the  fourth  example  neither  performed  very  well. 
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APPENDIX  A 

LISTINGS  OF  COMPUTER  PROGRAMS 


IDMCHKs  Control  routine 


EVEPS:,'  Subroutine  for  even-epsilon  algorithm 
EVRHO:  Subroutine  for  even-rho  algorithm 

TANDEM:  Merged  version  of  EVEPS  and  EVRHO 
SEQU:  Sequence  generators  for: 

(1)  s  -  ( 1+1 /n)n 

<>  ■  n 


(2)  s 


(3)  s 


II 

.•z<- 

k=»I 

n  •£" 


l)k+1/k 


(4)  sn  =  First  difference  of  logarithm  of  Turning  condition  number  for 

leading  n  x  n  segment  of  the  Hilbert  matrix.  Theory  from  Todd 

24  25 

(1954),  data  from  Savage  and  Lukacs  (1945). 

(5)  =  First  difference  of  logarithm  of  Turning  condition  number  for 

leading  n  x  n  segment  of  the  Hilbert  matrix.  Theory  from  Todd 

(1954),^  data  from  Fettis  and  Caslin  (1967). 
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i  .jw  ,»ty  *? 


PROGRAM  TOMCHK  73 /74  OPT«0  ROUNO«»-*/  TRACE  FTN  4.6*42 0 


1  PROGRAM  TOMCHKIXNPUT,  OUTPUT,  TAPE5  •  INPUT,  TAPE  6  »  OUTPUT) 

C 

d  EXERCISE  SUBROUTINE  TANOEM  FOR  CALCULATING  BOTH  EVEN  ORDERED 

C  EPSILONS  ANO  EVEN  OROEREO  RECIPROCAL  DIFFERENCES  FOR  VARIOUS 
9  C  INPUT  SEQUENCE S,  $N.  THE  LATTER  ARE  CEFINEO  BV  A  SEPARATE 

C  SUBROUTINE,  SEQU. 

C  f  ,  ■ 

C  THE  QUANTITIES  ENEN,  EOLO,  ANO  RNEH,  ROLD  ARE  FOR  USE  IN  TESTING 

C  FOR  CONVERGENCE,  IF  SO  OESIREO. 

10  C  I.  ■■ 

C  NS  IS  THE  NUMBER  OF  SKIPS  BEFORE  CALLING  TANDEM. 

C  NP  IS  THE  MAXIHUH  VALUE  OF  THE  ITERATION  INDEX,  NZ. 

C  NP  -  NS  =  MAXIMUM  VALUE  OF  INDEX,  IP.  ; 

C  " 

IS  C  IP  MUST  BE  SET  TO  ZERO  ON  INITIAL  CALL  TO  ASSURE  A  PROPER  START. 

ip  =  o  •  v  . 

NS  *  0 
NP  *  AO 
SHEW  =0.0 

20  ENEN  =  0.0 

RNEH  =  0.0 

10  HRITE<6,1010> 

1010  FORMAT <49H1 EVEN  ORDERED  EPSILONS  AND  RECIPROCAL  DIFFERENCES) 

11  HRITE (6, 1  Oil) 

25  1011  FORMAT (5  9H  SEQUENCE  SN  «  SUN ( i /N** 2) . . .. . <  (PI ) ♦* 2) /6  *  1.644934066 

IANS) 

12  NRITE<6, 1012) 

1012  FORMAT (39H0CVEN  OROER  *  2K  FOR  IP  =  2K*1  AND  2K*2//3M  IP,20X,5HIN P 
1UT , 1  AX , 7HEPSILON , 22X , 3HRH0/) 

30  15  DO  SO  NZ  =  1 ,NP 

EOLO  =  ENEM 
ROLO  =  RNEH 
SOLD  =  SNEH 

20  CALL  SEQU(NZ,SOLO,SNEM> 

35  IF<NZ  ,LE.  NS)  GO  TO  50 

SN  =  SNEH 

25  CALL  EVEPStSNjEP, IP,NS  ,IR) 

CALL  EVRHO(SN,RH, IP, NS , IR) 

ENEM  =  EP 

40  RNEH  =  RM 

IF (IR  .EQ.  1)  GO  TO  100 

30  HRITE<  6, 1030)  IP,SN,EP,RH 
1030  FORMAT  <1H  , 12, 3E25 . 14) 

5  0  CONTINUE 

45  100  HRITE (6,1100) 

1100  FORMAT (9HOFINISHEO) 

STOP 

END 


SUBROUTINE  EVEPS  >3/74  OP>80  ROUND8 ♦-*/  TRACE 


FJN  4. 6  ♦•*20 


SUBROUTINE  EVEPS <SN,EP , I P, NS ,  IR)  ..V  :  ,. 

EVEPS  COMPUTES  AN  ARRAY  OF  EVEN  OROEREO  EPSILONS,  Rll«>,  FOR  A 
GIVEN  SEQUENCE*  SN,  USING  THE  RECURSION  RELATION  GIVE N  BY  PETER  WYNN 
FOR  THE  EPSILON  ARRAY  AND  THE  ,PADE  TABLE. 

(NUN  NATH  9  (1966)  264-269) 

DIMENSION  Hl(SQ)  ,H2<50) ,Rtt50)  ,R2( 50)  ,Vl (50 >  ,V2C.50>  < 

IF  2*KHAX  IS  THE  HIGHEST  ORDEREO  RECIPROCAL  DIFFERENCE  DESIRED, 

THEN  THE  OIHENSION  OF  THESE  ARRAYS  MUST  BE  AT  LEAST  KMAX  ♦  1. 

THE  HIGHEST  ORDER  OF  EPSILON  CALCULATED  ANO  SENT  TO  OUTPUT  IS 
2K  FOR  IP  *  2K*1  ANO  2K*2. 

IP  MUST  BE  SET  TO  ZERO  ON  INITIAL  CALL  TO  ASSURE  A  PROPER  START. 

IR  =  0 

IR  IS  AN  ERROR  INDICATOR  WHICH  IS  SET  TO  1  IN  CASE  OF  PENDING 
DIVISION  BY  ZERO.  IT  MUST  BE  TESTEO  AND  PROPER  ACTION  TAKEN  IN 
THE  CALLING  PROGRAM. 


TOL  8  1.0E-14 

10  IP  8  IP  ♦  1 

IF (IP-2)  12,11,20 

11  V2(l)  8  1 .J /(SN  -  R2( 1 > > 

12  R2 (1)  8  SN 
H2  (1 )  8  0.0 
GO  TO  100 

21  JM AX  8  ( IP* 1) /2 
R1  (1)  8  SN 
HI (1 )  8  0.0 
00  26  J8 1 ,  JM  AX 
MS  8  Rl(J)  -  R2(J) 

IF (ABS (HS )  .LT.  TOL*  GO  TO  95 
21  VI  (J)  8  1.0/MS 
MS  8  H2t  J) 

IF ( J  . EO.  1)  MS  8  0.0 
23  HI ( J* 1 )  8  MS  ♦  VI ( J)  -  V2( J) 

IF  (ABS  (HI (J*l) )  .LT.  TOL)  GO  TO  95 
25  R1  (J  +  l)  8  R2(J)  ♦  1.0/HI ( J+l) 

3C  EP  8  Rl(JMAX) 

40  DO  43  J*1 , UMAX 
HZ  ( J )  8  HI  ( J) 

R2(J>  8  R1(J) 

43  V2  ( J)  8  Vl(J) 

GO  TO  100 

95  MRITE(6, 1095)  IP,J 

1095  FORMAT (26HODIVISOR  TOO  SMALL  AT  IP  8,I3,6H,  J  »,I2) 
IR  8  1 
100  RETURN 
END 


SUBROUTINE  EVRHO  73/74.  OPT*C  ROUND**-*/  TRACE  FTN  4.6*520 


10 


15 


20 


25 

t- 

30 


35 


40 


45 


50 


SUBROUTINE  EVRHO (SN.RH ,I>,NS, IR> 

£  .  «  "  .  .  .,  ;  •  .  ..  • 

C  EVRHO  COMPUTES  AN  ARRAY  OR  EVEN  0R0ERE0  RECIPROCAL  DIFFERENCES » 

C  RK.lt  FOR  A  GIVEN  SEQUENCE,  SN,  USING  A  RECURSION  RELATION  ANALOGOUS 
C  TO  THAT'  GIVEN  BY  PETER  MYNN  FOR  THE  EPSILON  ARRAY  ANO  THE  PACE  TABLE. 

C  (NUN  MATH  0  (1966)  264-269) 

C 

DIMENSION  Hl(50)  *H2(50)  ,R1(50)  »R2(50)  ,y.t<50)  ,V2(50) 

C  IF  i^KNAX  IS  THE  HIGHEST  OROEREO  RECIPROCAL  DIFFERENCE  DESIRED, 

C  THEN  THE  OIHENSION  OF  THESE  ARRAYS  MUST  BE  AT , LEAST  KM AX  ♦  1. 

C 

C  THE  HIGHEST  ORDER  RECIPROCAL  DIFFERENCE  CALCULATED  ANO  SENT  TO  OUTPUT  IS 
C  2K  FOR  IP  *  2K*1  ANO  2K*2. 

C 

C  IP  MUST  BE  SET  TO  ZERO  ON  INITIAL  CALL  TO  ASSURE  A  PROPER  START. 

C 

IR  *  C 

C  IR  IS  AN  ERROR  INDICATOR  WHICH  IS  SET  TO  1  IN  CASE  OF  FENDING 

C  DIVISION  BY  ZERO.  IT  MUST  BE  TESTEO  AND  PROPER  ACTION  TAKEN  IN 

C  THE  CALLING  PROGRAM. 

C 

TOL  «  l.Ct-14 
IF (IP-21  12,11,20 

11  V2 (1 )  =  1. 0/ (SN  -  R2 ( 1 ) ) 

12  R2 (1 )  -  SN 
H2  (1)  =  u.C 
GO  TO  10 0 

20  JM AX  =  (IP*l)/2 
Rl(l)  *  SN 

Hi  ( 1 1  *  C.G 
00  25  J=1,JMAX 
NS  *  R1(J)  -  R2(J> 

IF  (A6S(HS) '  .LT  •  TOD  GO  TO  95 

21  V1«J>  =  ( 2* J-l ) /MS 
MS  a  H2(J) 

IF (J  .EQ.  1)  MS  -  0.0 
23  H1(J*1)  =  MS  ♦  Vi(J)  -  V2(J> 

IF (ABS (H 1 ( J*l> )  .LT.  TOU  GO  TO  95 
25  Rt ( J*l>  *  R2 ( J )  ♦  (2*U)/H1(U*1> 

30  RH  a  R1 (UMAX) 

40  DO  43  J= 1 , JMAX 
H2  <  J)  *  H1IJI 
R2  ( J  >  *  R1(J) 

43  V2(J)  =  VI (J) 

GO  TO  100 

95  WRITE ( 6, 1C95)  IP,J 

1095  FORMAT (26HOOIVISOR  TOO  SMALL  AT  IP  =,I3,6H,  J  »,I2) 

IR  e  1 
100  RETURN 
END 
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SUBROUTI  HE  TANOEH  73/  74  OPTxO  ROUNOx,.*/  TRACE  FIN  4.B»420 


1  SUBROUTINE  TANOF  HI SN.EP, RHil I* ,  IR) 

C 

C  TANOEH  COMPUTES  TWO  ARRAYS  SIMULTANEOUSLY  FROM  A  GIVEN  INPUT  SEQUENCE 

C  SM,.  EVEN  ORtERfO  EPSILONS  IN  RE1  ARE  COMPUTED  USING  A  RECURSION 

S  C  RE  LAI  I  ON  GIVEN  PV  PETER  WYNN  FOR  THE  EPSILON  ARRA*  AND  PAOE  TABLE. 

C  (NUN  MATH  A  (19R6)  264-269) 

C  EVEN  OROEkEU  RES  IPROCAL  DIFFERENCES  IN  RR1  ARE  CONPUTEC  USING  AN 

C  ANALOGOUS  RECURSION  RELATION  DEVELOPED  BY  R.  P.  EDDY.  (UNPUBLISHED* 

■■■■■■■•  c  r 

to  DIMENSION  HE1  (50  > ,  HE?  ISO)  J  RE  1  (SO  1 ,  REKSQI,  VEK50I,  V E2(50) 

01  NENS10N  HR  1  ( SO  )  i  HR2  (SO  )  ,  RRU5  0  * ,  RR2(5Q>,  VRI  (SOI  ,  VK2(SQ) 

C  -  ■  ■  ••  •.  •  . 

C  RELATIONS  AMONG  PARAMETERS,.... 

C  LET  NO  BE  THE  DIMENSION  OF  THESE  ARRAYS.  AND  LET  2*KH»X  Bfc  THE 
IS  C  HIGHEST  ORDERED  EPSILON  AND/OR  RECIPROCAL  DIFFERENCE  UESIREQ. 

C  THEN  NO  .GE.  KNAX«1 

C  FURTHERMORE.  MAX  IP  x  NP-NS  =  2»NKAX»2  .LE.  2*N0 

C„  EXAMPLE...  NO  *  50,  2»KNAX  =  40,  MAX  IP  =  100 ,  NP  x  100,  NS  •=  0 

'c  '  '  .  .  . 

20  C  CONVERSELY,  FOR  A  GIVEN  IP,  TH.E  HIGHEST  ORDER  OF  EPS1LCN  AND/OR 

C  RECIPROCAL  0 1 FFERE  F'CE  CALCULATED  AND  SENT  TO  OUTPUT  VS  . 

'  C  2K  FOR  IP  =  2K«1  AND  2K»2, 

C 

C  IP  HUST  BE  SET  TO  ZERO  ON  INITIAL  CALL  TO  ASSURE  A  PRuPER  START, 

25  C 

IR  s  C 

C  IR  IS  AN  ERROR  INDICATOR  WHICH  IS  St  T  ll)t  t  IN  CASE  OF  FENDING 

C  DIVISION  BY  ZERO.  IT  MUST  BE  TESTED  AND  PROPER  ACTION  TAKEN  IN 

C  THE  CALLING  PROGRAH. 

30  C 

TOL  =  1.0E-14 
1C  IP  x  tp  .  i 

IP  (IP-2).  12,11,20 

11  VE 2 ( 1 )  *  1.0/  (SN  -  RE2<1>  I 

36  y  VR2(  1 )  x  VE2  ( 1 ) 

12  RE  2(1)  x  SN 

RR2I1I  x  SN 
HE?<1I  x  0.0 
HR 2 ( 1 1  x  0.0 

40  GO  TO  100 

?C  JNAX  X  (IPrl)/2 
RE  1  ( 1 )  x  SN 
RR 1  ( 1 )  x  SN 
HE  1  ( 1  )  x  0 . 0 

46  HR  1(1)  X  C.J 

00  25  J=1 , JMAX 

WE  x  PEl(J)  -  RE  2f  J ) 

WK  x  RRl(j)  -  RR  2 1 J ) 

IF  (ABS  (WE )  .LT.  Tol  .OR.  ABSSWK)  .LT.  TOL  I  GO  TO  96 
50  21  VE1U)  x  1.0/WE 

VRKJ)  x  <2*j-l)/WR 
WE  =  HE2(j) 

WR  x  mR2(j) 

IF  ( J  .NE.  II  GO  TO  23 
65  WE  x  0,0 

WR  X  0.0 

23  HE  1  ( J*  1)  =  WE  ♦  VE  1  <  J)  -  VE2IJ) 
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SUBROUTINE  TANDEM 


T3/7H  OP  T*  0  ROUNO*p-»/  TRACE 


FTN  A.6PA20 


HRHJ*1)  =  HR  ♦  VRKJ)  -  VR2  <  J) 

IF  (ABS  (HE  l(Jt  1))  .IT.  T0L  .OR.  ABS(HRlUtl)) 
RE1(J»1>  *  RE 2 (J t  ♦  1.C/HEKJUI 
25  RRHJ  +  l)  *  RR2JJ)  ♦  (2*J) /HR  1  ( j*l > 

30  EP  «  REl(JMAX) 

RH  a  RRl(JHAX) 

AO  00  A3  J*1 i JHA X 
HE2IJ)  s  HEi(J> 

HR2 ( J)  a  HR1 (J) 

RE2CJI  «  REKJ) 

RR2( J)  a  RRl ( j) 

VE2I J)  «  V£1 ( J) 

A3  VR  ?  ( J )  a  VKUJ) 

CO  TO  100 

95  HR ITE  1 6> 1 095 )  IP,J 

1095  FORMAT (26H001 VISOR  TOO  SNAU  AT  IP  s,I3,6H» 
IR  =  1 
IOC  RETURN 
END 


.LT.  TOU  CO  TO  95 


J.  >. 12) 


Bliliii'liiiliiBiL 


o  o 


SUBROUTINE  SEQU 


7  3/74 


OP  T* 0  POUNO  *♦- * /  TRACE 


FTN  4.6+4ZJ 


01/14/77 


1  SUBROUTINE  SE QU< N7 , SOLO, S NEK , SN) 

PARTIAL  SUNS  OF  SUM < < ( -1 > •* ( N -II » /N> . LN(Z>  *  0.6931471E055B945 

SN  IS  SNITCH  TO  EFFECT  ALTERNATING  SIGN  IN  SUMMATION. 

J.C  SN  a  -SN 

11  SNEM  =  SOLO  ♦  SN/NZ 
10C  RETURN 


5 


o  o  o 


SUBROUTINE  SEQU 


23/2 4  OPT«Q  ROUNO**-*/  TRACE 


FtN  4.6*420 


t 


5  10 

ig  %> 


SUBROUTINE  SEQU(N2»SOLD»SNEM> 

PARTIAL  SUNS  SN  ■  SUMU/N»*2>  FROM  l  TO  M2. 
LIMIT  IS  tP!**2W6  »  1.64493  40666  46... 

SNEN  *  SOLO  *  l.O/NZ**2 
RETURN 


TPHirV 3r':'  vV’Pi^'  n'  "*'  ■"  •- -r{’^ T. r" 


SUBROUTINE  SEQU  7 J/74  0»T*C  ROUND**-*/  TRACE  FTN  4.6*420 


1  SUBROUTINE  SE  QUt  NT .SOLD* SNEH) 

C  ASVHPIOTIC  SPECTRAL  CONDITION  NUMBER  OF  LEADING  SEGMENTS  OF  THE 

C  HILBERT  MATRIX . FXP(3.525*NI 

DIMENSION  EKLHlOh  EVLN<10> 

5  VL1<  l)  *  1.0 

EVll<  21  *  2. 0/3.1  (SORT  (1  3.0))  /6.  0 
EVL  1  <  3)  *  1.4083189271237  E0 
EVL1  <  4>  «  1 . 500  21 42  8  0  05  9  2  EO 
EVLH  5)  *  1 • 567  050691 0982  E0 
10  EVL 1 (  6)  s  1.6188998589243  EQ 

EVL1 (  7)  =  1.6608853399269  E3 

EVL1 1  8)  ■  1.6959389969219  EO 

EVL1 (  9)  s  1.7258826609018  EO 

EVLt  ( 1 0)  ^  1.751  9196702652  EO 

15  EVLN  (  1)  =  1.0 

EVLNI  2)  =  2. 0/3.0  -  ( SORT  (1  3  .OH  /6.  0 

EVLN  <  3)  =  2.6873403557735  E-Q3 

EVLN<  4)  =  9 . 67C  23 C40 2258 7  E-05 

EVLH<  5)  =  3.2879287721719  E-06 

20  EVLNI  6)  s  1.0827994845656  E-07 

EVLN (  7)  s  3.4938986059912  E-09 

EVLN (  8)  *  1.1115389663724  E-10 

EVLN <  9)  =  3.4936764029115  E-12 

EVLN ( 1 0)  =  1.0931538193797  E-13 

25  1C  A  =  <EVLl(NZ>*EVLN<N2-l))/»EVLN(N2>*EVLi(N2-in 

11  SNEH  =  AL  OG ( A ) 

100  RETURN 
ENO 


45 


APPENDIX  B 
SAMPLE  OUTPUT 

Output  pages  show  the  performance  of  EVEPS  and  EVRHO  on  each  of  the  sequences 
given  In  Appendix  A.  The  first  three  sequences  are  well  known,  as  are  their  limits: 

Cl)  e  »  2.71828  18284  59045 

(2)  ln2  -  0.69314  71805  59945 

(3)  ir2/6  -  1.64493  40668  48226 

24 

(4)  Following  Todd  one  can  show  that,  for  large  n,  the  Turning  condition 
number  of  leading  segments  of  the  Hilbert  matrix  is  asymptotically 

CT  -  (l/4ir2/2)  EXP(nB) 

where 

B  -  2 In  ((✓Y+l)/(/2-D)  -  3.525496 
This  B  is  the  limit  of  the  sequence  s^. 

(5)  It  is  qonjecturable  that,  asymptotically,  the  spectral  condition  number  is 

1*  ..  i  '■ 

the  same  as  the  Turning  condition  number  for  leading  segments  of  the  Hilbert 
matrix.  This  seems  to  be  borne  out  by  numerical  evidence  from  the  even-rho 
algorithm. 


47 


oaoaaiina 
o  □  a  o  o  U  lj  :3 

111  III  III  III  III  111  Ui  111 

ajua'fNJ  p«« 
aa'fw^MU'lTi^ 
C5»HS.f>»PO»t5(OvO 
O  iD  Ni  IMP  N  IB  4 

cs^u5t*>i-3»-ir0(\i 

Of~i«N.»***NN 

a»W«INaM<Ji 

oi-(K>crcT''tjiaN» 

O  ifi  O'  N  O  iC  lO  4* 

O  JIPNNO'OO 
O'  LP  «  IP  *  IP  'O  vO 


0  00000000 
a  o  o  a  a  13  15  a  o 
+  »  +  #  +  »  +  +  + 
UJUJUJUJUJUJUJIUtjJ 
N  n  a  ^  Ni  if\  irs  n  O'* 
HNWO'IPiS'H'04 
M  >5  K)  o  ip  N  «  ^  |0 
M  ifl  •»  «(p  Q  MO  a 
NWIMPiOIPtOStO 
v£  OHfllp^lPrt4 
i«IOkrl»  NNN 
4  »l(Pa 

4«0'tf'vjDN.M0'<v» 

OOlP'JSl^r-IIOvOlT'JO 
Ifi  rl  (O  S  4  O'  ITi  O' 
(T'lPtOPjtS«rOO)K. 


FINISHEO 


Ei/EN  ORDEREO  EPSILONS  AND  RECIPROCAL  DIFFERENCES 

SEQUENCE  SN  =  SUM<1/N**2> . (<PI>**2»/6  «  1.644934066848^4 

EVEN  OROER  *  2K  FOR  IP  »  2K*t  ANO  2K*2 


IP 


INPUT 


EPSILON 


RHO 


1  .  lOOOQOOOOOOOOOE+Ol 

2  • 125000000000006*01 

5  . 13611illiniiiE*oi 

**  *142  36111111111 E*  Cl 

5  • 146361111111116*01 

6  • 14913880888889E*01 

7  . 1511797Q521542E*01 

•  . 15274220521542E*01 

9  .15397 6773116 65 E*01 

10  • 154976773116656*01 

11  • 155803219397656*01 

12  . 15649 7663042096*01 

13  . 157Qd937981842E*01 

14  • 157599503900056*01 

19  .158044028344506*01 

16  .158434653344506*01 

17  . 15878C67413574E*Cl 

18  .159089316081056*01 

19  . 15936632439130E+01 

20  .159616324391306*01 

21  .159843081760926*01 

22  .160049693331166*01 

23  .  160238729247996*01 

24  .  160412340359106*01 

25  .160572340359106*01 

26  .160726269353186*01 

27  .160857443564436*01 

28  . 160984994584846*01 

29  .161103900649056*01 

30  . 16121 5011760 16E*C 1 

31  .161319070032796*01 

32  . 161416726282796*01 

33  . 16150855364735E+01 

3 V/  .161595058837666*01 

35  .161676691490726*01 

36  .161753851984556*01 

37  ■ 161B268980Q354E+Q1 

38  .161896150081106*01 

39  • 161961 8963 uC 69 6*01 

40  . 162024396300696*01 

finished 


0. 

0. 

.145000000000006*01 
.150396825396836*01 
. 155161 74402250£+6i 
. 15 71 76736854 75E+ 01 
.1590 305 413615  6E*G1 
.15999841551 50  IE *01 
.160908690629166*01 
.161447429518606*01 
. 161960991390246*01 
.162291529061566*01 
.16260 947416967E*01 
« 162826822009076*01 
.16  3  0  3  72  5  5  23  8  746  *01 
. 16 318 771 93290 OE* 01 
. 16333452531 30 6E *01 
. 163441739948236*01 
. 163553392 650 28E *01 
• 16362259559 1406*01 
.163686943326696*01 
.163575972445046*01 
•  1637971546892 4E* 01 
• 16 3744044 78456E* 01 
.1637701190714  76*01 
.163853748696396*01 
.16306790 62590 5E* 01 
.163859734257846*01 
.16 386565 82639 4E *01 
. 163854346c 0 99 6E* 01 
• 1 6 3955863 158 2 4E* 01 
. 163937343647646*01 
. 16 3948 75 100 90 7E *01 
.164060083105376*01 
. 16410161 8 620 69E* 01 
.164064070784146* 01 
.  16401 1470  89  31  OE'- 01 
. 16 402 9 54 5 7881 3E* 01 
. 16404572 4 888 9 4E* 01 
.  16397048229620E*fll  , 


.0. 

0. 

. 165000000900006*01 
.164682539682546*01 
. 1644894  8948949E*  01 
.1644922586521 IE* 01 
.164493437612376*01 
.164493414537346*01 
. 164493406439676*01 
. 16 44934 06621 93 E* 01 
. 1644934067341 9E* 01 
. 164 4934066562 7E*  31 
. 1644934 06691 35E* 01 
. 164493406676546*01 
. 16 4 49340 66999 4E* 01 
•  164493406680096*01 
. 16449340 66841 8 E*  01 
. 164493406600 14E+01 
.  164 49340 6684 7 OE* 01 
V 16  44934  0668398 E* 01 
. 16449340668470  E*  31 
. 164493406684286*01 
. 164493406684666*01 
. 164493406690046*01 
.164493406604766*01 
. 16 4 4934 06661 8 5E* 01 
.164493406684636*01 
. 1644934 066S457E* 01 
. 16 4 493406684 64 E* 01 
. 1644934 0 668429 E*  31 
. 1 6 449340 66 8462 E*  01 
.164493406684626*01 
• 16449340668462E*0i 
. 164493406685196*01 
. 16449340668462E*01 
.164493406684656*01 
. 164493406684926*01 
.164493406684646+01 
. 1644934 06684 63E+  01 
. 164493406684656*01 
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